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Summary 


Deconvolution  and  estimation  of  transfer  function  phase  and 
coefficients  for  nonGaussian  linear  processes 


>  1 

i 


n 


NonGaussian  linear  processes  are  considered.  It  is  shown  that 
the  phase  of  the  transfer  function  can  be  estimated  under  broad  condi¬ 
tions.  This  is  not  true  of  Gaussian  linear  processes  and  in  this  sense 
Gaussian  linear  processes  are  atypical.  The  asymptotic  behavior  of  a 
phase  estimate  is  determined.  The  phase  estimates  make  use  of 
bispectral  estimates.  These  ideas  are  applied  to  a  problem  of  decon¬ 
volution  which  is  effective  even  when  the  transfer  function  is  not 
minimum  phase.  A  number  of  computational  illustrations  are  given. 
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1.  Introduction.  Assume  that  the  random  variables 


v^.t  =  .  .  .  ,  -1, 0,  1, . .  .  are  independent  and  identically  distributed  with 

2 

mean  zero,  Ev  =  0,  and  variance  one  Ev  =  1.  Let  [a.  }  be  a  sequence 
t  t  j 

of  real  constants  with 


I 


a .  <  ® 
J 


Consider  the  linear  process  generated  by  { a. ^  }  and  {v^} 

CO 

■*\L  v-j  • 


(1) 


Let  a(z)  =  £  a  .  2?  be  the  z-transform  corresponding  to  the  process  {x 

j  J  * 


Then 


a(e 


•*S-S 

i  * 


is  called  the  frequency  response  function  or  transfer  function.  We  are 
concerned  with  the  estimation  of  a  (e  on  the  basis  of  observations 

only  on  the  process  {xfc}  . 

The  spectral  density  of  {x^  }  is 

«»  *  jjf  |o(e'U)  |  2  . 

In  the  Gaussian  case  (when  {x^}  is  normally  distributed)  the  full  proba¬ 
bility  structure  of  {x^}  is  determined  by  f ( A )  or  equivalently  by  the 
modulus  of  a(e  |  a(e  *  ^ |  .  The  phase  information  in  a(e  is 


i 
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not  identifiable  in  the  Gaussian  case. 


If  a(z)  is  a  rational  function 


a(z) 


A(z) 

B(z) 


-with  A(z),  B(z)  polynomials 


q 

A(z)  =  S  ak  *k  • 

k  =  0 


P 

■  X  \  *k  • 

k  =  o  k 


o  . 


=  i  , 


the  process  [x^  }  is  a  finite  parameter  autoregressive  moving  average 
process,  that  is, 

p  q 

V«-j =  |0v..k  •  <» 

If  {xt}  is  a  Gaussian  process  satisfying  (2),  then  any  root  z.  4  0  of  A(z)  or 
B(z)  can  be  replaced  by  its  conjugated  inverse  z’1  without  changing  the  proba¬ 
bility  structure  of  {x^}  ,  This  follows  since  jel^-z0|  =  | zQ| "2 |e**- z‘ 1  |  . 

If  all  the  roots  are  distinct  there  are  2P+C*  ways  of  specifying  the  roots 
without  changing  the  structure  of  fx^  .  To  ensure  unique  determination  of 
the  coefficients  a^  and  b^  of  (2)  (since  there  is  a  different  specification  of 
these  coefficients  corresponding  to  each  of  the  2p+q  root  specifications)  in  the 
Gaussian  case,  it  is  the  custom  to  assume  that  all  the  roots  of  A(z)  and  B(z) 
are  outside  the  unit  circle  |z|  <  1  in  the  complex  plane. 


ftir- 


T 


However,  for  a  nonGaussian  stationary  process  satisfying  (2)  (in  j 

p+q 

which  case  the  independent  vt's  are  nonGaussian)  the  different  2 

specifications  of  roots  mentioned  above  generally  correspond  to  different  j 

probability  structures  and  different  processes.  As  a  simple  example, 
consider  the  moving  average 


xt  *  6vt  -  5vt-l  tVt-2 


with  the  roots  o'  A(z)  2  and  3,  and  the  moving  average 

yt  =  3vt  -  7Vl  +  ^t-2 

having  a  polynomial  A(z)  with  roots  1/2  and  3.  Both  {x^  }  and  f  } 
have  the  same  spectral  density  but  if  the  independent  random  sequence 
{ v^  }  is  exponentially  distributed,  the  marginal  distributions  of  the  {xt} 
and  {yt  }  sequences  are  different.  In  the  problem  of  deconvolution  where 
one  wishes  to  recover  the  process  {v^}  (assumed  nonGaussian  which  is 
most  often  the  case  in  applications)  in  some  sense,  the  proper  specification 
of  roots  (which  are  inside  and  which  are  outside)  becomes  crucial  (see 
Rosenblatt  (1974)).  There  is  a  discussion  concerning  the  distribution  of 
roots  as  related  to  prediction  problems  in  Rosenblatt  (1980). 

There  are  results  on  the  estimation  of  the  coefficients  a.  and  b 

J  k 

of  (2)  (corresponding  to  roots  assumed  outside  the  unit  circle)  in  Box  and 
Jenkins  (1976).  In  the  Gaussian  case  these  are  essentially  equivalent 
asymptotically  to  maximum  likelihood  procedures.  In  the  nonGaussian  case 
the  computations  are  carried  out  as  if  the  process  were  Gaussian.  One  has 
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a  least  squares  but  not  a  maximum  likelihood  solution  in  the  nonGaussian 
case.  The  coefficients  estimated  are  those  corresponding  to  roots  outside 


the  unit  circle  even  though  the  actual  structure  of  the  process  may  not  be 
one  with  all  the  roots  outside  the  unit  circle.  Thus  one  will  typically  not  be 
able  to  resolve  the  actual  structure  using  these  procedures  in  the  non¬ 
Gaussian  case.  Of  course,  if  one  knows  the  actual  nonGaussian  distribu¬ 
tion  of  the  v  ,  one  can  use  the  maximum  likelihood  estimate  or  an  asymp¬ 
totically  equivalent  procedure  to  estimate  the  coefficients  in  (2)  even  if  the 
roots  z.  are  not  all  outside  the  unit  circle  (see  Bhasawa,  Feigen  and 
Heyde  (197  5)).  Higher  order  spectral  methods  discussed  in  the  next  section 
do  not  require  this  knowledge.  Our  discussion  follows  that  of  Rosenblatt 
(1980). 


2.  Higher  order  spectral  method 

Assume  that  {x^}  is  a  linear  process  (see  (1)),  with  the  independent 
random  variables  f  v^  }  nonGaussian  and  having  all  moments  finite.  Actu¬ 
ally  we  only  require  that  some  cumulant  of  order  k  >  2  be  nonzero. 
Also  let 


and  assume 


2  1 3 1  I  a. |  <  - 
j  1 


a(e_i  *  )  4  o 

for  all  X.  We  will  see  that  a(e  is  essentially  identifiable  (as  con¬ 


trasted  with  the  Gaussian  case)  if  one  only  observes  the  process  [x^  }  . 


wrw 


Since  the  v^'s  are  assumed  nonGaussian  with  all  moments  finite, 
there  must  be  a  cumulant  of  v  ,  y  i  0  of  smallest  subscript  k  >  2.  The 

i  K 

k  order  cumulant  spectral  density  of  the  process  }  is  given  by 


VS . xk-i)s:rVi.  2  .  cum(’vvi. . xt«.  > 


<2b>  j, . V-i 


(-!?!«•>.) 


Vl  “(«  1)"-a(e  k'‘)  o( 


•Uk-l\  I 
e  /ale 


Then 


h(  X)  =  arg 


i.  -iX.  a 

a(e  >U 


a  (l) 
a(l)i 


) k  ^-1  * 

taOii)  Yk=,2">  bk(0 . 0)/ t«0)32 


h(  X  j)  +  •  •  •  +h(  X  k  +  j) 


7 log  [' 


xl\  ad)  )k  -1 


\\ail)\]  ^  b^Ul . Xk-1) 


Cf(  Xj)  •  •  •  f(  X  k_j)  £(  X  x  +  •  •  •  +^k.1)}‘1/2 


^»V**J*’  .-•4  - 


■  'Mr*#*:***'  i 


since 


Further 


h(-X)  =  -h(  X)  . 


h'(0)-h'(X)  =  lim  — ~ r  f  h(X)  +  <k- 2)  h(A)  -  h(  X  +  (k-2)A)] 
A  -»  0 


Now 


(5) 


h(  X)  =  J  {  h/(u)  -  h/(0)  }du  +  c  X  =  h  (  X)  +  c  X 
0 


where 


In  particular 


c  =  1/(0)  . 


h(TT)  =  h^(  tt)  +  err  . 


Since  the  a.'s  are  real  we  must  have 
J 


h(  tt)  =  k  TT 


for  some  integer  k.  Set 


Then 


so  that 


h1(TT)/TT  =  6  . 


h(n)  =  krr  =  (6  +c)tt 


c  =  k  -  5  . 


The  integer  k  cannot  be  determined  without  further  assumptions  since  it 
corresponds  to  reindexing  or  subscripting  the  v^'s.  The  sign  of  a(l)  is 
also  intrinsically  undecidable  since  one  can  multiply  a^'s  and  v^'s  by  (-1) 
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without  changing  the  observed  process  {x^}  •  Thus,  under  the  conditions 


specified  above  for  a  nonGaussian  linear  process  {xt }  ,  a  (e  1  *  )  is 
identifiable  up  to  the  integer  k  and  the  sign  on  the  basis  of  observations 
on  {x^}  only  and  is  given  by 

a(e_i^)  =  |  2TTf(  X)|  ^  2  exp  { ih( X)  } 


f\ 

h(X)  =  J  {h  /(u)  -  h'(0)}du  +  cX 


h  (tt) 

=  hj(X)  -  -rs—  ^  - 

Notice  that  h^(X)  can  actually  be  computed. 

3.  Phase  estimation  and  convergence  of  estimates 

There  are  many  discussions  concerned  with  the  estimation  of  the 
second  order  spectral  density  f(  X)  (see  Anderson  (1971)  or  Jenkins  and 
Watts  (1968)).  We  will  concentrate  on  the  estimation  of  h(  X).  For 
simplicity  of  discussion  we  will  assume  that  the  third  order  cumulant  y  ^ 
of  v  is  nonzero.  The  program  in  the  higher  order  case  can  be  carried 
out  in  a  similar  manner.  Equation  (5)  becomes 


h'(0)-h'(X)=  lim  fh(X)  +  h(A)  -  h(  X  +A)} 
A  -»  0  ^ 


when  k  =  3.  For  (4)  we  find  that  up  to  a  sign 

h(  X)  +  h(  A)  -  h(  X  +  A)  =  arg{b3(X,  A)}  . 

From  this  point  on  we  will  drop  the  subscript  and  understand  that  we  are 


mu-  -***•»> 
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dealing  with  the  bispectral  density  b(X,n).  Let  b(X.u)  be  an  estimate 

n 

of  b(  X  ,  (i)  based  on  a  sample  of  size  n.  Then  an  estimate  of 


9(  X  ,  u)  =  arg  b(  X  ,  M) 


can  be  given  by 


9  (  X  ,  u)  =  arctan(Im  b(X,u)/Re  b(  X  ,  u)) 
n  n  n 


We  note  that  for  a  complex  number 


z  =  x  +  iy  =  re 


i9 


with  r  =  |  z \  and  9  =  arctan(y/x)  a  principal  value  determination,  one  has 


99  _  x  99 

9y  2  ’  9x 
r 


9  29  _  2xy  9^9  _  2xy 

a  2  '  4  ’  2  "  4 

9x  r  9y  r 


and 


Therefore 


9^9  _  _1_  2x2  _  1  2y  2 

9x  9y  2  4  2  2 

r  r  r  r 


(X.  u)  -  9(X,  u)  =  -  -m  b(  X  -  U l  {Re  b(X.  u)  -  Re  b(  X  ,  u)  } 

lb(X.U>i2 

+  — -  b^  -  {Im  b(X,u)  -  Imb(X.U)} 
|b(X,  [X)\ 


+  o  (  b(  X  ,  u  )  -  b(  X  ,  n))  . 
P  n 


(7) 


10 


as  an  appropriate  estimate  of  h ^ (X ).  Assume  that  sixth  moments  are  finite 
and  that  bispectral  estimates  ^b(X  ,  | a)  of  the  type  obtained  by  weighted 
averages  of  3rt*  order  periodogram  values  (see  Brillinger  and  Rosenblatt 
(1967))  are  employed.  It  has  been  shown  in  the  paper  just  cited  that  if  the 
bispectral  density  is  continuously  differentiable  up  to  second  order  and  one 
has  a  symmetric  bandlimited  weight  function  with  bandwidth  A,  that  then 
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^  "  b(X  ,u)  +  vD  )  b(X  ,  u)  w(u,  v)  dudvA  2  +  o(A  2) 

K* 


and 


rr2(nb(X,u))  ~  f<  X  +  H>  /w2(u.v)dudv  . 


An 


(9) 


(10) 


if  A  n->® 


as 


n  -»  ®  ,  A(n)  -►  0.  Here  and  represent  partial  derivative 


with  respect  to  X  and  ^  respectively.  Further  estimates  b(X,u), 

n 

nb(X  >U  )t  0  ^  p,  <  X  ,  0<(j.  5  X  ,  are  asymptotically  uncorrelated  as  n 
if  (X  »|i)  and  (X  #u  )  are  distinct.  Using  (7)  and  (8)  we  can  write 

Hn(  X)  -  hx(  X)  =  Rn(  X)  +  op(Hn(  X)  -  hj<  X  )) 

and  show  by  employing  (9)  that 
.X 


ER  (X) 


n 


l 


Im  b(u,  0)  j-AD  2  Re  0)  +  2B  D  D  Re  b(u,  0) 
"  u  v 


[0  )  b(u,  0)  | 


2  ‘  u 


+  CDv  Re  b(u,  0)]  duA 


(ID 


i 


+  I  — e  bfu,  y  [AD  2  Im  b(u,  0)  +  2BD  D  Im  b(u,  0) 

C  U  U  V 


b(u,  0)  | 


+  CD^  Im  b(u,  0)]  du  A 

+  o(A) 

where  A,  B,  and  C  are  the  second  moments  of  the  weight  function  w  of 
the  bispectral  estimates 
2 


A  =  Ju  w(u,  v)  du  dv  ,  B  =  J  uvw(u,  v)dudv 
C  =  J  v2w(u,v)dudv  . 
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Further,  by  using  (10)  it  follows  that 


cov(R  (X),  R  (u))  - 
n  n 


m  f 


min(X,|i) 


f f ^(u)/ |b(u,  0) |  Z]  du  w^(u,v)dudv.  (12) 


Assuming  the  existence  of  all  moments  and  (11),  one  can  show  that  H  (X) 

n 

is  asymptotically  normal  with  mean  h^(X)  and  variance  given  by  (4  2).  The 
mean  square  error  of  R^(X)  is  of  order 

A  2  C3 

c>  rr 

A  n 

-2/5  -1/5 

and  the  optimal  rate  of  convergence  is  n  when  A(n)  ~  n 

Assume  that  the  bispectral  density  is  continuously  differentiable 
up  to  third  order.  Further  let  the  weight  function  of  the  estimate  be  band- 
limited  with  first  and  second  order  moments  zero.  Such  a  weight  function 
cannot  be  nonnegative  everywhere.  Then 


E  b(X  ,u)  -  b(X  ,  u)  =  0(A  ). 
n 


The  mean  square  error  of  R^(X)  is  of  the  order 


a4  +  C  2 

ciA  +~r  • 

A  n 

The  optimal  rate  of  convergence  is  then  n  4^  with  A(n)~  n’*^  . 

Generally  we  will  estimate  h^(X)  and  hence  h(X)  for  a  whole  range 
of  X  values.  The  sign  of  b(0,  0)  may  not  be  positive.  We  estimate  it  by 
noting  the  real  part  of  ^bfA.A).  If  it  is  negative  we  multiply  all  ^bfjA.A) 
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r 


with  a  minus  sign.  The  estimate  H^(X)  is  then  given  by 

k-  1 

H  (X)  =  -  ^  argf-  b(jA,  A) }  . 

n  n 

J  =  1 

4.  Computations  using  spectral  methods 

We  remark  on  the  computational  aspect  of  phase  estimation  of 
a(e  and  give  a  few  illustrative  examples  to  indicate  its  effectiveness. 

Given  a  sample  [x^]  °f  size  n  =  kN,  we  center  and  normalize  it  so 
that  it  has  mean  zero  and  variance  one.  Break  up  the  sample  into  k  dis¬ 
joint  subsections  of  equal  length  N  so  that  the  variance  of  the  bispectral 
estimate  from  each  section  is  not  too  large.  Then  choose  a  grid  of  points 

X  .  =  jA  in  (0,  2tt),  j  =  1 . M,  A  =  2ttL/N.  Though  the  symmetry  condi- 

3 

tion  h(X)  =  -h(-X)  implies  that  one  need  only  deal  with  X  in  (0,rr),  there 

may  be  some  advantage  in  considering  X  €  (0,  2tt).  We  will  comment  on 

this  point  later  on.  Form  the  bispectral  estimate  ^bfjA.A)  of  the  type 

discussed  above  with  a  weight  function  of  bandwidth  A  from  each  subsection. 

Average  the  estimates  from  the  different  subsections  so  as  to  arrive  at  a 

final  estimate  b(jA,A).  A  detailed  discussion  of  this  kind  of  algorithm  can 
n 

be  found  in  Helland  and  Lii  (1981).  Compute  9  (j)  =  arg  {  b(jA,A)}  +  2kn 

n  n 

where  the  integer  k  is  chosen  to  ensure  continuity  of  H^(f  A)  =  H^(X^)  = 
i  .  i 

-  Zj  9  (j),  l  =2 . M+l  (neighboring  values  are  as  close  to  each  other 

as  possible).  Since  the  upper  index  is  l  -1  we  start  with  1=2.  Since 

h( 0)  =  0  one  sets  H  (0)=0  and  estimates  H  (A)  =  H  (X  )  by  an  interpolation 

n  n  n  i 

between  0  and  H  (X ,),  X  =  2A.  H  (tt)  is  also  computed  by  an 
n  &  2  n 
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interpolation  procedure.  This  amounts  to  a  complete  procedure  for 


estimating  h(  X). 


Since 


i  r  , 

ak  =  2 n  J  tt< 


-i\.  ikX 
e  )  e  dX 


and  estimate  of  is  given  by 


1  f  -i 

’j's!  a(e 


-iX.  ikX  ,. 
)e  dX 


M  +  1 

2TT(M+ij  .^0  '/2TTf 


_  (  ,  H  (TT)  J 

(  X  .)  exp  i  H  (X  .)  -  X  .  +  kX.  H 

n  j  n'  j  tt  j 


and  this  computation  can  be  carried  out  by  using  the  fast  Fourier  transform. 

The  a^'s  are  real  numbers  and  so  the  a^'8  may  or  may  not  be 
real.  If  the  symmetric  property  of  f(  X)  and  h(  X)  is  used  and  the  integra¬ 
tion  is  carried  out  from  -tt  to  tt  almost  real  a.  's  will  be  obtained.  The 

k 

imaginary  part  of  the  ct^'s  will  only  be  the  size  of  rounding  errors.  In 

practice  there  is  no  indication  of  how  good  or  bad  the  estimates  are  apart 

r  i  M+l 

from  asymptotic  results.  In  actual  practice  [jA.L_Q  may  not  be  sym¬ 
metric  about  tt  .  If  the  estimates  H^(Xj)  are  reasonably  good  the  estimated 
a/s  from  (13)  should  still  be  almost  real.  The  size  of  the  imaginary  part 

reflects  the  level  of  noise.  When  the  estimates  H  (X  .)  are  not  good  the 

n  J 

imaginary  part  of  the  a^'8  becomes  comparable  to  (or  larger  than)  its  real 
part.  This  can  serve  as  a  direct  indication  of  the  quality  of  the  estimation. 


wm 


If  the  linear  process  is  one-sided  with  a  finite  number  of  parameters 


one  has  a  moving  average  of  order  q 


+  0  . 


q  4  q  ^  ; 

We  could  estimate  a(z)  =  £.  q  cl.  zj  by  a(z)  =  a. .  zJ  .  In  deconvolu- 

J  J  3 

tion  we  try  to  recover  the  process  [v  }  ,  v  =  x.  (B  is  the  backward 

t  t  a(B)  t 

shift  operator  so  that  x  =  x  .)  by  computing  the  approximation 

t  t-j 

-  1 

v.  =  x.  .  If  all  the  roots  of  a(z)  (and  a(z))  are  outside  the  unit  circle 
t  a(B)  t 

(the  frequency  function  is  minimum  delay)  then  a  *(z)  has  a  one-sided 

as  . 

expansion  £j_Q  a  /  B^  .  In  the  computation,  the  series  is  truncated  after  a 
certain  number  of  terms.  If  some  of  the  roots  of  a(z)  have  modulus  less 
than  one  we  can  still  expand  &  *(B)  with  a  Laurent  series  expansion.  Once 
the  roots  of  a(z)  are  computed,  one  can  easily  get  the  Laurent  series 
expansion  of  a  *(B)  by  partial  fractions  as  described  in  Rosenblatt  (1974) 
or  Henrici  (1974). 


Another  way  to  find  the  inverse  weights  in  deconvolution  is  to  use 
a  least  squares  criterion  as  described  in  Wiggins  (1978).  Another  general 
method  of  deconvolution  will  be  mentioned  in  the  section  on  computation. 


5.  Other  possible  computational  methods 


We  briefly  discuss  two  other  possible  methods  of  estimating  the 
coefficients  of  a  nonCaussian  moving  average  process  of  order  q  . 


a 

t.  =  X  a .  v  . 
1  j  =  0  3  4-3 


As  noted  earlier,  second  order  moments  will  not  allow  us  to  determine  the 
location  of  the  roots  of 


a(z)  =  ^  a.z3  . 
j  =  0  3 


Higher  order  moments  will  be  used  in  the  first  method  which  makes  use  of 

_  3 

a  least  squares  procedure.  Assume  Evt  s  0,  Evt  =  y  ^  0.  Consider 


ci=EVt.k 


Xaj  aj+k  '  k  =  -q,-q+l . . 


Estimate  by 


sk4  £  v 


and  solve  the  extremal  problem 


7  X(^h^)  • 


There  are  q+2  unknowns  aQ,  .  . .  ,  a  and  y  in  (17).  Due  to  the 


homogeneity  of  the  a/s  we  have  to  normalize  the  problem  appropriately; 
all  the  a.'s  can  be  multiplied  by  a  constant  c  0  and  y  can  be  divided  by 
c^  without  changing  (17).  There  are  a  number  of  ways  of  carrying  out  such 
a  normalization.  One  could  set  Ev^  =  y  =  1.  Alternatively  a^  =  1  could 
be  the  normalization  condition.  Some  comments  on  the  asymptotic  distribu¬ 
tion  of  the  c^'s  are  given  in  Appendix  2. 

The  second  method  is  a  searching  procedure.  One  uses  a  typical 


second  order  method  to  estimate  the  roots  of  a(z),  r  ^  ,  j  =  1,  . .  .  ,q,  assuming 
all  the  roots  have  modulus  greater  than  one.  An  accurate  estimate  of  the 
distribution  of  roots  is  obtained  by  taking  the  conjugated  inverse  of  an 
appropriate  number  of  the  r/s.  Suppose  all  of  the  r/s  are  real  and  dis¬ 
tinct.  Then  there  are  2^  possible  sets  of  roots  that  give  the  same  second 


order  structure.  Each  of  these  sets  yields  a  distinct  set  of  a/s  which  in 

turn  lead  to  a  distinct  set  of  the  c,  's.  Choose  the  set  of  a.'s  which  deter- 

k  J 


mine  the  set  of  c^'s  minimizing 


among  all  the  possible  sets  of  {  cL  }  .  If  some  of  the  roots  r^  are  complex, 
the  inverse  complex  conjugates  are  taken  in  pairs.  If  there  are  multiple 
roots,  the  solution  of  roots  in  terms  of  coefficients  is  unstable.  Some  com¬ 
ments  on  this  question  are  made  in  Appendix  1.  The  initial  set  of  coeffi¬ 
cient  estimates  corresponding  to  roots  all  outside  the  unit  circle  can  be 
obtained  by  the  method  described  in  Box  and  Jenkins  (1976).  Alternatively, 
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F 


one  could  try  to  obtain  the  roots  directly  by  solving  for  the  roots  of  the 
polynomial 

p(z)  =  zq  g(z ) 

where 

g(z)  =  a(z)  a(z"  S 


with 


We  estimate  [3.  by 


8.  zj 

J 


3i  .  i  =  E  x.  x 


lit  *t  +  |j| 

a<  a*+ui 


n*  ijl 

2,  xtxt+m 


The  roots  of  p(z)  with  modulus  greater  than  one  are  the  initial  set  of  roots. 


6.  Examples 


We  will  consider  a  few  simple  examples  generated  by  Monte  Carlo 
simulation  to  illustrate  the  computation  and  to  give  a  qualitative  feeling  of 
the  effectiveness  of  the  theory.  Details  and  possible  "fine  tuning"  of  the 
computational  method  will  be  considered  elsewhere. 

We  generate  xt  =  vt  +  vt  j  +  a2  Vt  2  '  *  =  *  *  •  »  ^40  where 

v  =  v'  -  1  and  v^'s  are  independent  exponentially  distributed  random 
deviates  with  mean  one  obtained  from  the  GGZEN  subroutine  in  the  inter¬ 
national  Mathematical  and  Statistical  Library  (IMSL).  Then 


E  v  =  0  ,  »  i  •  2v.  =1  » 

t  t 


Evt  =  2  . 
r  1  640 

We  partition  {.x^i^^  into  five  lections,  each  of  which  has  128  points. 
Compute  the  bispectrum  estimate  b(*L(jA,  A)  ,  j  =  1,  13;  i  =  1, . .  .  ,  5  ,  by 
the  algorithm  described  in  Lii  and  Helland  (1981).  Here  we  set 
A  =  y|g  =  0.  442.  Our  final  bispectrum  estimate  is 


5 

b640«A-A>.*  I  S  b<,a<JA.A>  • 

i  =  1 


Compute 


9  (j)  =  arg{b  (jA.A)} 


n 


n 


=  arctan(Im  b  (jA,A)/Reb  (jA , A )) 


n 


n 


by  taking  the  principal  value  as  well  as 

j-  1 


-H  (jA)  =  -H  (X.)  =  7  9  (i)  ,  j  =  2,  ...»  14  . 

n  n  J  i=!  n 
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Let 


H  (0)  s  H  (XJ  =  0  , 

n  n  0 


H  (A)  =  H  (X.)  =  H  (X  )  . 
n  n  l  d  n  d 


H  (n)  =  i[H  (7A>  +  H  (8A )] 
~  d  n 


n 


n 


=  (3.  093)  +  H  (3.  534)] 

2  n 


n 


and 


c  =  -H  (tt)/tt  = 
n 


-  6  . 


Recall  that  c  is  an  estimate  of 

..  h(A) 
c  =  lim  — — 

A  -i  0  A 


h'(0) 


up  to  an  integer.  We  will  use  c  =  h(A)/A  instead  of  h'fO)  to  compare 

with  c  in  the  following  examples.  From  formula  (9)  we  compute 

A  standard  smoothed  periodogram  with  uniform  weights  and  bandwidth  A 

was  used  to  compute  f  (X)  as  an  estimate  of  the  spectrum  f (X )  of  fx  ]  (f  (0)  is 

n  t  n 

obtained  by  a  linear  extrapolation).  These  examples  are  as  follows: 


Model:  =  vt  +  a  ^  vt  j  +  a2vt  2  ^our  c'*ses  specified  given  below: 


Case 


1 

2 

3 


Coefficients  Roots 


a  0 

al 

a2 

ri 

r2 

1.0 

-0.833 

0.  167 

2.0 

3.  0 

1.0 

-2.  333 

0.  667 

0.  5 

3.0 

1.0 

-3.  50 

1.  50 

2.0 

0.  333 

1.  0 

-5,0 

6.  0 

0.  5 

0.  333 

4 


Case  1.  c  r  1.29,  c  -  1.605 


X 

Length 
!a(e  **)! 

Est. 

Length 

Argument 

h(X) 

Argument 

by  Sum 

-H  (X) 
n 

Argument 

at  J 

9  (X) 
n 

Adjusted 

H  (X)  +  c\ 
n 

0. 

0. 3333 

0.  181 

0. 

0. 

0. 

0. 

0.442 

0.4194 

0.  399 

0.5732 

0. 3125 

0. 

C. 3966 

0.884 

0.6509 

0.617 

0.8309 

0.6249 

0.6  249 

0.  7933 

1.  325 

0.  9776 

1.  050 

0.8428 

1. 3195 

0.6945 

0.8078 

1. 767 

1. 3393 

1.452 

0. 7180 

0. 9188 

0. 5994 

0.  9176 

2.209 

1.6685 

1.663 

0.5199 

2. 7462 

0.8273 

0.  7994 

2.651 

1. 9032 

1.  715 

0.2830 

3.  5893 

0.  8431 

0.6654 

3.093 

1.9990 

2.241 

0.0286 

4.4842 

0.8949 

0.4795 

3.5  34 

1.9376 

1.  740 

-0.2274 

5.6009 

1. 1166 

0. 0720 

3.  976 

1. 7307 

1.702 

-0.4697 

6.2542 

0.6534 

0.  1277 

4.418 

1.4176 

1.455 

-0.6788 

7. 1482 

0.8939 

-0.0571 

4.860 

110573 

1.  174 

-0.8240 

7.9014 

0.7532 

-0. 1012 

5.  301 

0.7172 

0.655 

-0.8502 

8.  1520 

0. 2506 

0. 3573 

5.  743 

0.4599 

0.445 

-0.6584 

9.  3127 

1. 1607 

-0.  0943 

6.  185 

0. 3377 

0.235 

-0. 1461 

9.4044 

0. 0917 

0.  52  30 

P> 

o 

=  0.  9593 

• 

o 

1 

tl 

H 

<d 

5816  ,  a3  = 

0. 1158 

Here  and  from  this  point  on  all 


are  adjusted  by  sign  and  index  shift. 
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Case  2 


c  =  -1.  388  ,  c  =  0.4094 


■  X 

Length 

|a(e 

Est. 

Length 

Argument 
h(  X) 

Argument 

by  Sum 

-H  (X) 
n 

Argument 

at  J 

9  (X) 
n 

Adjusted 

H  (>  )  +  cX 
n 

Hn(X, 
+(e  -  1)X 

0.6667 

0.473 

0. 0000 

0. 

0. 

0. 

0. 

0.8389 

0.  894 

-0.6125 

-0.2713 

0. 

0.0904 

-0.  352 

B 

1.  3018 

1.  315 

-1.0828 

-0.5425 

-0. 5425 

0.  1808 

-0.  703 

1.  325 

1. 9553 

1.888 

-1.4915 

-0.4851 

0.  0574 

-0.0575 

-  1. 382 

1.767 

2.6785 

2.  748 

-1.8895 

-0.6904 

-0.  2053 

-0.0331 

-1.80 

3.  3369 

3.  348 

-2.2893 

-0.6481 

0.  0423 

-0.2563 

-2.465 

2.651 

3.8064 

3.  520 

-2,6920 

-1. 0097 

-0.  3615 

-0. 0756 

-2. 727 

3.9980 

-  3.  0966 

-1.0884 

-0.  0787 

-0.  1778 

-3.  27 

3.  534 

3.8752 

3.  951 

-3.  5014 

-1.4842 

-0. 3958 

0.  0371 

-3.497 

3.  976 

3.4614 

3.297 

-3.  9047 

-1.6797 

-0. 1956 

0.0518 

-  3.  924 

4.418 

2.8351 

2.802 

-4. 3051 

-1.5538 

0.  1259 

-0.2550 

-4.693 

4.860 

2.1146 

2.  168 

-4. 7031 

-1.8428 

-0.  2890 

-0.  1469 

-5.007 

5.  301 

1.4344 

1.429 

-5. 1070 

-1.8785 

-0.  0357 

-0.2921 

-5.593 

5.  743 

0.9199 

0.  936 

-5. 5558 

-2. 3478 

-0.  46  93 

-0.  0037 

-5.  747 

6.  185 

0.6755 

0.443 

-6. 1 366 

-2.6185 

-0.  27  06 

0.  0861 

-6.099 

a 

a 

0 

=  0. 7164 

.  o1  = 

-2.175,  a3 

=  0.  7605 

As  an  indication  of  discretization  error,  notice  that  if  we  use  the  exact  h(X)  and 
|n  (c  *  )|  instead  of  estimated  ones,  we  get 

aQ  =  0.9136,  a1  =  -2.247,  a3  =  0.  5977 
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Case 

3.  c  =  -0. 

613,  c  ; 

=  0.6363 

Argument 
at  J 

e  (>) 

n 

Adjusted 

H  (>)  +  cX 
n 

H  (X) 

n'  ' 

+  (c  -  1)X 

X 

Length 

|a(e‘iX)| 

Est. 

Length 

Argument 

h(X) 

Argument 

by  Sum 

-H  (X) 
n 

0. 

1. 0000 

0.  055 

0 

0 

0. 

0 

0. 

0.442 

1.2583 

1.280 

-0.2711 

0. 0062 

0. 

0.2749 

-0.  167 

0.884 

1. 9527 

2.  505 

-0.6843 

0. 0125 

0.  0125 

0. 5498 

-0,  334 

1.  325 

2. 9330 

2.  922 

-1.  1592 

0. 3266 

0.  3141 

0.5168 

-0.808 

1.767 

4.0178 

4.  386 

-1.6448 

0.5819 

0. 2553 

0.5426 

-1.224 

2.209 

5. 0054 

5.272 

-2. 1286 

1.0098 

0.  4280 

0.  3958 

-1.813 

2.651 

5.  7097 

6.  732 

-2.6094 

1. 1805 

0.  1707 

0.  5062 

-2. 145 

3.093 

5.  9971 

5.440 

-3.0884 

1.8244 

0.6439 

0.  1435 

-2.  95 

3.  534 

5.8129 

6.670 

-3.5672 

2. 1738 

0. 3494 

0. 0752 

-3.459 

3.  976 

5. 1922 

5.610 

-4. 0475 

2.5575 

0. 3836 

-0. 0273 

-4. 003 

4.418 

4.2527 

4.653 

-4.5306 

2.8902 

0.  3328 

-0. 0790 

-4.497 

4.860 

3.  1720 

3.  040 

-5.0162 

2.9963 

0. 1061 

0.0961 

-4. 764 

5.  301 

2. 1516 

2.630 

-5.4959 

3.  3573 

0.  3610 

0.  0161 

-5.215 

5.  743 

1.  3799 

1.456 

-5.9306 

3. 3021 

-0.0552 

0.  3525 

-5. 391 

6.  185 

1. 0132 

0.282 

-6.2334 

3. 3388 

0.  0367 

0. 5970 

-5. 588 

a 

ao 

=  0. 7561 

,  ax  =  -3. 

334  ,  a2 

=  1.  778 

9 


a2  =  1.  778 


Case. 4 

.  c  =  -  3. 

2  9,  c  - 

1.  18 

Argument 

by  J 

9  (X) 
n 

Adjusted 

H  (X)  +  aX 

n 

H  (X) 
n 

+(c  -  2)X 

X 

Lerr*h 

1  a(e  1?l)  ! 

Est. 

Length 

Argument 

h(X) 

Argument 

by  Sum 

-H  (X) 
n 

0. 

2. 0000 

0.  969 

0 

0 

0. 

0 

0. 

0.442 

2.  5167 

2.403 

-1.4567 

0. 0542 

0. 

-0.  5757 

-1.4597 

0.884 

3.  905  3 

3.  837 

-2.5981 

0. 1085 

0. 1085 

-1. 1514 

-2. 9194 

1.  325 

5.8659 

5.  832 

-3.4935 

-0.  5056 

-0. 6140 

-1. 0589 

-3.  7089 

1.  767 

8.  0356 

6.  112 

-4.2523 

-1. 0588 

-0. 5533 

-1. 0270 

-4. 561 

2.209 

10. 0109 

9.601 

-4. 9377 

-1. 7068 

-0. 6480 

-0.  9006 

-5. 3186 

2.651 

11.4194 

10. 311 

-5. 5844 

-2. 5605 

-0.  8537 

-0.5683 

-5.  8903 

3.  093 

11. 9942 

10. 165 

-6.2136 

-3.2706 

-0.7101 

-0. 3797 

-6. 5657 

3.  5  34 

11.6258 

10.822 

-6.8412 

-4. 1459 

-0.  8753 

-0.0259 

-7.  0939 

3.976 

10. 3843 

9.  952 

-7.4825 

-4.  9293 

-0.7834 

0.2360 

-7.  716 

4.418 

8.  5055 

6.  902 

-8. 1569 

-5.6188 

-0. 6895 

0.4040 

-8.432 

4.860 

6. 3440 

5.  925 

-8.8953 

-6.0445 

-0. 4257 

0.  3083 

-9.412 

5.  301 

4. 3033 

4.  360 

-9.7527 

-6.6180 

-0.  5735 

0.  3603 

-10.242 

5.743 

2. 7597 

2.560 

-10.8280 

-6.7895 

-0.  1716 

0.0104 

-11.495 

6.  185 

2.0265 

0.  760 

-12.2239 

-6.8955 

-0. 1060 

-0.4051 

-12. 775 

a 

ao 

=  0.  9603 

,  a1=-3. 

904,  &2  = 

4.  966 

When  we 

increase  the  number  of  parameters  to  four. 

we  observe 

qualita- 

tively  the  same  type  of  result  as  in  the  three  parameter  case.  But  generally 
speaking,  the  instability  is  increased.  We  give  one  example:  the  model  is 

xt =  vt  -  4-25  Vt-1  +4-75  Vz  •  °-  938  Vi 

2  1 

with  roots  —  ,  ^ ^  ,  4.0,  we  get: 
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C  =  -3.6887  ,  c  =  -1.825 


X 

Length 

|afe-iX)| 

Est. 

Length 

Argument 

h(X) 

Argument 

by  Sum 

-H  (X) 
n 

Argument 
by  J 

e  (X) 

n 

Adjusted 

H  (X )  +  cX 
n 

H„(M 

+(c  -  2)X 

0. 

0.5625 

0. 

0. 

0. 

0. 

0. 

0. 

0.442 

0.  9470 

0.  925 

-1.6304 

-0.  3495 

0. 

-0.4569 

-1.  341 

0.884 

2.  0233 

2.202 

-2.6631 

-0.6990 

-0. 6990 

-0.  9138 

-2.682 

1.  325 

3.7839 

3.439 

-3.4612 

-1. 1051 

-0.  4062 

1 

H- < 

• 

UJ 

M 

o 

-3.  964 

1.  767 

6. 0495 

6.  768 

-4. 1782 

-2. 1441 

-1. 0389 

-1.0814 

-4.615 

2.209 

8.  3719 

8.  149 

-4.8647 

-3.2201 

-1.0760 

-0.8118 

-5.23 

2.651 

10.  1662 

11.574 

-5.5391 

-4. 0550 

-0.8349 

-0.  7832 

-6. 085 

3.093 

10.  9296 

10. 548 

-6.2089 

-5.0367 

-0. 9818 

-0.6079 

-6.  794 

3.  5  34 

10.4383 

11. 975 

-6.8783 

-6.4316 

-1. 3949 

-0.0193 

-7.  087 

3.  976 

8.8369 

8.413 

-7.5511 

-7.2046 

-0.7730 

-0.0527 

-8. 005 

4.418 

6.5805 

7.  820 

-8.2339 

-8.  1175 

-0.9129 

0.0538 

-8.  782 

4.860 

4.2557 

3.  946 

-8.9415 

-8.8852 

-0.7677 

0.0152 

-9.  7U5 

5.  301 

2.  3572 

2.  518 

-9.  7138 

-9.  7556 

-0.8704 

0.  0792 

-10.  >2 

5.743 

1.  1270 

1.  160 

-10.6702 

-10.0662 

-0. 3106 

-0.4166 

-11.  9 

6.  185 

0.5829 

0. 

-12. 1458 

-10.  7554 

-0.6893 

-0.  5337 

-12.  0 

A 

ao 

=  1. 106, 

=  -3. 

572,  a2  = 

4.293,  a3 

=  -0.  9762 

26 


These  simple  examples  indicate  that  one  can  estimate  the  unknown 

coefficients  reasonably  well  and  one  is  able  to  discriminate  different 

models  even  though  they  have  the  same  spectral  structure.  The  question 

of  determining  how  many  roots  are  inside  of  a  unit  circle  can  be  answered 

very  reliably  using  this  method.  One  simply  takes  the  absolute  value  of  an 

estimate  H  (X)  +  c(X)  of  h(X)  near  X  =  2tt  and  divides  it  by  2rr.  round- 

ing  the  result  to  its  nearest  integer.  This  integer  is  the  winding  number 
-iX 

given  by  a(e  )  which  gives  the  number  of  zeros  of  a(z)  inside  of  a  unit 
circle.  One  can  see  that  this  is  clearly  the  case  in  all  these  examples. 

Graphs  l  through  5  compare  the  theoretical  h(X)  with  the  estimate 
obtained  by  using  our  techniques  in  the  five  successive  cases  considered. 
Graphs  6  through  10  are  concerned  with  deconvolution.  We  consider  the 
moving  average 

x.  =  v.  -  2.  333  v,  ,  +  0.  667  v 
t  t  t-1  t-2 


with  the  v  's  independent  exponential  variables  with  mean  one.  The  v  's 

»  X 

are  generated  as  pseudo-random  variates.  F(x)  is  the  exponential  distribu¬ 
tion  function  with  mean  one.  Graph  6  is  a  plot  of  the  sample  distribution 

2 

function  of  the  generated  v^'s  and  of  F(x).  Let  a(B)  =  1  -  2.  333B  +0.  667  B  . 
v  is  obtained  by  truncating  the  exact  deconvolution  formula  as  applied  to 
the  generated  time  series 


A 


9 


1 

a(B)  xt  * 


(18) 
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Graph  7  has  a  plot  of  the  sample  distribution  function  of  the  v^'s  with  F(x). 
The  estimated  model  (using  our  bispectral  techniques)  is 

x.  =  0.7164  v  -  2.  17  5  v.  +  0.7605  v.  ,  . 

t  t  t  *  I  t*  M 

2  si 

Let  a(B)  =  0.7164  -  2.  17  5  B  -  0.7605  B  .  Let  v.  =— -  x,  with  the  same 

a(B)  t 

truncation  as  that  employed  in  (18).  Graph  8  is  a  plot  of  the 
sample  distribution  function  of  the  with  F(x).  In  these  three  graphs 
a  one  sample  95%  confidence  band  using  the  Kolmogorov- Smirnov  statistic 
is  indicated.  Graph  9  has  plots  of  the  sample  distribution  functions  of  the 
v^'s  and  v^'s  respectively.  Graph  10  plots  the  sample  distribution  functions 

A 

of  the  vt's  and  the  v^'s.  Two  sample  95%  confidence  bands  are  given  on 
these  two  graphs. 

A  general  way  to  find  the  deconvolution  weights  can  be  described 
as  follows.  We  have  an  estimate. 

a(e  )  =  2.  Tr'f  (X )  exp  {i  (H  (X)  +  c  X)  )  . 

n  n 

To  find  a  (e  ** )  we  compute  b(e  **  )  =  exp  { -0n(  a  (e  **))  }  .  Then 

TT 

,  -iX  ikX  . 
b(e  )  e  dX 

give  the  coefficients  of  series  expansion  of  a  ^ (e  **)  which  are  the  weights 
desired  in  deconvolution.  This  method  is  general.  We  do  not  require 
knowledge  of  the  order  of  the  moving  average  process  and  there  is  no 
need  to  compute  the  roots  from  the  estimated  coefficients  and  how  the 
roots  are  distributed  is  irrelevant. 
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Appendix  1 


In  this  appendix  we  consider  the  relationship  between  the  coeffi¬ 
cients  of  a  polynomial  and  the  roots  of  the  polynomial,  at  least  locally. 


The  polynomial  is 


X  a  zj  =  [{  (*-»  ) 

j=0  3  j=l  3 


where  the  roots  are  z.  ,  j  =  1, . .  . ,  p,  and  the  coefficients  a.  ,  j  =  0,  1 . . 

with  a  =1.  It  is  well  known  that 
P 


1  =  -2  zj  ’ 


ap"2  jTk  ZjZk 


i  1  z .  z  z 

p-3  jA*/  1  k  ' 


a0  -  (-1}  zj  •  •  •  zp 


Let  us  consider  the  relationship  between  the  differentials  of  the  coefficients 

a.  and  the  differentials  of  the  roots  z,  .  Now 
J  k 


9a  _  v 

=  Z,  z-  =  _a  i  ‘  z#  » 

azi  fri  1  p-1  1 


1 


1 


1 


1 


V  = 


U  is  a  triangular  nonsingular  matrix  and  V  is  the  Vander  Monde  matrix. 
V  is  nonsingular  as  long  as  the  roots  z.  are  distinct. 
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Appendix  2 


The  object  is  to  remark  on  some  aspects  of  the  asymptotic  behavior 
of  estimates  of  the  c^'s  in  the  context  of  a  general  linear  process.  Let 
x^  be  a  nonGaussian  linear  process  (1)  with 

Ev  =  0  E  \  “  =  1  Ev^=y 

t  *  t  ’  t  y3 

cum(v^)  =  y4  ,  cum(vS  = 

Set 


yt=Xt  ‘ 


For  convenience  we  introduce 


g  -  co  v(x  ,  y  )  , 
u  t  t-u  ’ 

r  =  co  v(x  ,  x.  )  , 
u  t  ’  t-u 

h  =  cov(y  ,  y  )  . 
u  7t  t-u 


Consider  the  estimates 


N 

K =  b  2  xt  yt +, 


of  Ex  vx,  .  Then 
t  t+  a 


N 


c°v(g  ,  g.)  =  \  ^  eov(x  y  x  y 

a  hT  t,T  =  1  t  t+a  T  T 


+  b 
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where 


COv(xtyt+a*  Xr  yT+b) 


rt-T  Ct-T  +a-b  +  -b  ®t-T+a 


+  cum(xt,  xT,  yT+b)cum(yt+a) 
+  cum(xt,  xT,  yt+a)cum(yt+b) 

+  cum(xt,  yt+a.  xT,  yT+b)  . 


Now 


r  =  /  a.  a. 

u  4*  j  j-u 


V  2 

g  =  y  ^  ?  o  .  a. 

Bu  3  4  j  j-u 


h  =V,^ctZa2  +  2  ( S  a.  a.  \ 
u  4  ^  j  j-u  \4-  j  j-u / 


Further 


cum(yt+a)  =  ^  a2  , 

j 

cum(xt,  xt,  yT+b)  =  Y4  Xajaj2+baj+(t.T)=k<b.t-T)  . 
cum(xt,  xT,  yt+a)  =  Y4Sa.a2+a  oj+(T_t)  =k(a,x-t)J  , 


cum(xt,  yt+a,  xT,  Yj  +  b>  =  Y&  £  au  a 

u 


2  2 
u+(T-t)  u+a  u+b+(T-t 


+  4yA  (  7  a  a 

4\1r  u 


u+(x-t)  ttu+a  au+b+  (t 


au  ttu+  (T-t)+b-aj 
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f 


+  2y3  (S  au  au+(T-t)  au+a  )  (2au  au+(t-T)+a-b) 

+  2y  3  au  au+(t-T)  au+bj  au  au+(T-t)  ib  zj 

=  k(a,b,t-T)  . 

It  is  then  clear  that 

lim  N  cov(ga.  i^)  =  J  r  c  +  2  B„+a 

N  ->  00  s  8 

+  rQ  ^  [k(b,  s)  +  k(a,  s)  3  +  ^  k(a,b,  s)  . 
s  s 

Under  the  assumption  that  E  I  a ^  ]  <  ®  ,  with  a  truncation  argument  like 
that  employed  in  Anderson  (1970),  one  can  show  that 

/N(ga-  ga) 

are  asymptotically  normally  distributed  with  covariance  structure  given 
by  the  preceding  formula. 
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